

## #######################################################################################################
## 
## Simulation Replication R command file for:
## Graeme Blair and Kosuke Imai, "Statistical Analysis of List Experiments," Political Analysis, In Press.
##
## NB: The main replication file for the tables and figures is in this archive as PA_Replication.R
##
## Created 27 October 2011
## 
## Contact Graeme Blair <gblair@princeton.edu> with questions
## 
## #######################################################################################################

## ####
## Simulations of the bias, consistency, and coverage of alternative list experiment estimators

  ## #####
  ## Generate the simulations

rm(list = ls())
library(MASS)
library(gamlss.dist)
source("Corstange2009_listit.R")

## Create simulated data

library(list)

logit <- function(x) exp(x)/(1+exp(x))
inv.logit <- function(x) return(log(x)-log(1-x))

sim.num <- 0

sim.object <- list()

# set up simulation

## set up simulation
sims <- 10
ssize <- 500
alpha <- 0.1
	
for(y.prop in list(c(1/2, 1/2, 1/2, exp(.5)/(1+exp(.5))), c(1/2, 1/4, 3/4, exp(.5)/(1+exp(.5))), c(1/5, 2/5, 3/5, 4/5, exp(.5)/(1+exp(.5))), c(1/6, 3/6, 4/6, 4/6, exp(.5)/(1+exp(.5))) ) ){

	sim.num <- sim.num + 1
	
	print(paste("=== starting sim", sim.num))
	
	# set up true values
	
	# last covariate/intercept is treatment
	y.int <- rep(0, length(y.prop))
	
	y.cov <- log((y.prop)/(1-(y.prop)))*2
	
	par.control <- c(y.int[1], y.cov[1])
	for(k in 2:(length(y.int)-1)) par.control <- c(par.control, y.int[k], y.cov[k])
	
	par.treat <- c(y.int[length(y.int)], y.cov[length(y.int)])
			
	J <- length(y.int)-1
	
	par.control.list <- list()
	for(j in 1:J) par.control.list[[j]] <- par.control[(j*2-1):(j*2)]
	
	first.mle <- FALSE
	
	est.treat.ml <- se.treat.ml <- est.treat.nls <- se.treat.nls  <- est.treat.listit <- se.treat.listit  <-  cov <- matrix(NA, ncol = length(par.treat), nrow = sims)
	  
	est.control.ml <- se.control.ml <- est.control.nls <- se.control.nls <- est.control.listit <- se.control.listit <- matrix(NA, ncol = length(par.control), nrow = sims)
	
	for (i in 1:sims) {
		
		print(paste("starting sim",i))
	
		## data generation
		dat <- data.frame(x=runif(ssize))
		x <- model.matrix(~ x, data = dat)
	
		y <- list()
		for(j in 1:J){
			y[[j]] <- rbinom(n = ssize, 1, logit(x %*% par.control.list[[j]]))
		}
		y[[J+1]] <- rbinom(n = ssize, 1, logit(x %*% par.treat))
		
		y <- do.call(cbind, y)
		
		y.control <- y[,1:J]
		y.treat <- apply(y,1,sum)
		treat <- sample(rep(c(0,1), ssize/2), size = ssize, replace = FALSE)
		
		for(j in 1:ncol(y.control)) y.control[treat == 1,j] <- NA 
		y.treat[treat==0] <- NA
		
		rownames(x) <- NULL
					
		simdata <- data.frame(x=x, y.control, y.treat, treat)
		names(simdata) <- c("(Intercept)","x",paste("y.control", 1:ncol(y.control),sep=""),"y.treat", "treat")
						
		print("starting ml")
		
		sim.ml <- ictreg(as.formula(paste("cbind(", paste(paste("y.control", 1:ncol(y.control),sep=""),collapse=", "), ", y.treat) ~ x")), data = simdata, method = "ml", first.mle = first.mle)
		
		est.treat.ml[i, ] <- sim.ml$par.treat
		se.treat.ml[i, ] <- sim.ml$se.treat
		
		est.control.ml[i, ] <- sim.ml$par.control
		se.control.ml[i, ] <- sim.ml$se.control
		
		print("starting nls")
		
		sim.nls <- ictreg(as.formula(paste("cbind(", paste(paste("y.control", 1:ncol(y.control),sep=""),collapse=", "), ", y.treat) ~ x")), data = simdata, method = "nls", first.mle = first.mle)
						
		est.treat.nls[i, ] <- sim.nls$par.treat
		se.treat.nls[i, ] <- sim.nls$se.treat
		
		est.control.nls[i, ] <- sim.nls$par.control
		se.control.nls[i, ] <- sim.nls$se.control
		
		print("starting listit")
		
		x.noint <- x[,-1, drop = FALSE]
		
		mod <- listit(Y = cbind(y.treat, y.control), X = x.noint)
				
		est.treat.listit[i, ] <- mod$b.star
		se.treat.listit[i, ] <- mod$se
		
		est.control.listit[i, ] <- mod$b.control
		se.control.listit[i, ] <- mod$se.control
			  
	}
	
	## confidence intervals
	
	critical <- abs(qnorm(alpha/2))
	
	# coverage for ml
	
	lb.treat.ml <- est.treat.ml - se.treat.ml*critical
	ub.treat.ml <- est.treat.ml + se.treat.ml*critical
	
	cov.treat.ml <- matrix(NA, ncol = 2, nrow = sims)
	for(i in 1:sims){
		for(j in 1:length(par.treat)){
			cov.treat.ml[i,j] <- ((lb.treat.ml[i,j] < par.treat[j]) & (par.treat[j] < ub.treat.ml[i,j]))
		}
	}
	
	
	lb.control.ml <- est.control.ml - se.control.ml*critical
	ub.control.ml <- est.control.ml + se.control.ml*critical
		
	cov.control.ml <- matrix(NA, ncol = length(par.control), nrow = sims)
	for(i in 1:sims){
		for(j in 1:length(par.control)){
			cov.control.ml[i,j] <- ((lb.control.ml[i,j] < par.control[j]) & (par.control[j] < ub.control.ml[i,j]))
		}
	}
	
	# coverage for nls
	
	lb.treat.nls <- est.treat.nls - se.treat.nls*critical
	ub.treat.nls <- est.treat.nls + se.treat.nls*critical
	
	cov.treat.nls <- matrix(NA, ncol = 2, nrow = sims)
	for(i in 1:sims){
		for(j in 1:length(par.treat)){
			cov.treat.nls[i,j] <- ((lb.treat.nls[i,j] < par.treat[j]) & (par.treat[j] < ub.treat.nls[i,j]))
		}
	}
	
	
	lb.control.nls <- est.control.nls - se.control.nls*critical
	ub.control.nls <- est.control.nls + se.control.nls*critical
		
	cov.control.nls <- matrix(NA, ncol = length(par.control), nrow = sims)
	for(i in 1:sims){
		for(j in 1:length(par.control)){
			cov.control.nls[i,j] <- ((lb.control.nls[i,j] < par.control[j]) & (par.control[j] < ub.control.nls[i,j]))
		}
	}
	
	# coverage for listit
	
	lb.treat.listit <- est.treat.listit - se.treat.listit*critical
	ub.treat.listit <- est.treat.listit + se.treat.listit*critical
	
	cov.treat.listit <- matrix(NA, ncol = 2, nrow = sims)
	for(i in 1:sims){
		for(j in 1:length(par.treat)){
			cov.treat.listit[i,j] <- ((lb.treat.listit[i,j] < par.treat[j]) & (par.treat[j] < ub.treat.listit[i,j]))
		}
	}
	
	
	lb.control.listit <- est.control.listit - se.control.listit*critical
	ub.control.listit <- est.control.listit + se.control.listit*critical
		
	cov.control.listit <- matrix(NA, ncol = length(par.control), nrow = sims)
	for(i in 1:sims){
		for(j in 1:length(par.control)){
			cov.control.listit[i,j] <- ((lb.control.listit[i,j] < par.control[j]) & (par.control[j] < ub.control.listit[i,j]))
		}
	}
	
	sim.object[[sim.num]] <- list(est.control.ml=est.control.ml,  se.control.ml=se.control.ml, cov.control.ml=cov.control.ml, ub.control.ml=ub.control.ml, lb.control.ml=lb.control.ml, est.treat.ml=est.treat.ml, se.treat.ml=se.treat.ml, cov.treat.ml=cov.treat.ml, ub.treat.ml=ub.treat.ml, lb.treat.ml=lb.treat.ml, est.control.nls=est.control.nls, se.control.nls=se.control.nls, cov.control.nls=cov.control.nls, ub.control.nls= ub.control.nls, lb.control.nls=lb.control.nls, est.treat.nls=est.treat.nls, se.treat.nls=se.treat.nls, cov.treat.nls=cov.treat.nls, ub.treat.nls=ub.treat.nls, lb.treat.nls=lb.treat.nls, est.control.listit=est.control.listit, se.control.listit=se.control.listit, cov.control.listit=cov.control.listit, ub.control.listit=ub.control.listit, lb.control.listit=lb.control.listit, est.treat.listit=est.treat.listit, se.treat.listit=se.treat.listit, cov.treat.listit=cov.treat.listit, ub.treat.listit=ub.treat.listit, lb.treat.listit=lb.treat.listit, par.control=par.control, par.treat=par.treat)
		
}

  ## ####
  ## Figure 4 - "Monte Carlo Evaluation of the Three Estimators"
  ## pg 22

load("ModifiedN500.RData")
load("ModifiedN1000.RData")
load("ModifiedN1500.RData")

sims <- 10000

# SIMULATION 1

## bias

nls1.500.bias <- apply(nls1.500$est.treat[1:sims,], 2, mean) - c(0,1)
nls1.1000.bias <- apply(nls1.1000$est.treat[1:sims,], 2, mean) - c(0,1)
nls1.1500.bias <- apply(nls1.1500$est.treat[1:sims,], 2, mean) - c(0,1)

mle1.500.bias <- apply(mle1.500$est.treat[1:sims,], 2, mean) - c(0,1)
mle1.1000.bias <- apply(mle1.1000$est.treat[1:sims,], 2, mean) - c(0,1)
mle1.1500.bias <- apply(mle1.1500$est.treat[1:sims,], 2, mean) - c(0,1)

listit1.500.bias <- apply(listit1.500$est.treat[1:sims,], 2, mean) - c(0,1)
listit1.1000.bias <- apply(listit1.1000$est.treat[1:sims,], 2, mean) - c(0,1)
listit1.1500.bias <- apply(listit1.1500$est.treat[1:sims,], 2, mean) - c(0,1)

## rmse
nls1.500.rmse <- sqrt(nls1.500.bias^2 + apply(nls1.500$est.treat[1:sims,], 2, var))
nls1.1000.rmse <- sqrt(nls1.1000.bias^2 + apply(nls1.1000$est.treat[1:sims,], 2, var))
nls1.1500.rmse <- sqrt(nls1.1500.bias^2 + apply(nls1.1500$est.treat[1:sims,], 2, var))

mle1.500.rmse <- sqrt(mle1.500.bias^2 + apply(mle1.500$est.treat[1:sims,], 2, var))
mle1.1000.rmse <- sqrt(mle1.1000.bias^2 + apply(mle1.1000$est.treat[1:sims,], 2, var))
mle1.1500.rmse <- sqrt(mle1.1500.bias^2 + apply(mle1.1500$est.treat[1:sims,], 2, var))

listit1.500.rmse <- sqrt(listit1.500.bias^2 + apply(listit1.500$est.treat[1:sims,], 2, var))
listit1.1000.rmse <- sqrt(listit1.1000.bias^2 + apply(listit1.1000$est.treat[1:sims,], 2, var))
listit1.1500.rmse <- sqrt(listit1.1500.bias^2 + apply(listit1.1500$est.treat[1:sims,], 2, var))


## coverage

nls1.500.cov <- apply(nls1.500$cov.treat[1:sims,], 2, mean)
nls1.1000.cov <- apply(nls1.1000$cov.treat[1:sims,], 2, mean)
nls1.1500.cov <- apply(nls1.1500$cov.treat[1:sims,], 2, mean)

mle1.500.cov <- apply(mle1.500$cov.treat[1:sims,], 2, mean)
mle1.1000.cov <- apply(mle1.1000$cov.treat[1:sims,], 2, mean)
mle1.1500.cov <- apply(mle1.1500$cov.treat[1:sims,], 2, mean)

listit1.500.cov <- apply(na.omit(listit1.500$cov.treat[1:sims,]), 2, mean)
listit1.1000.cov <- apply(na.omit(listit1.1000$cov.treat[1:sims,]), 2, mean)
listit1.1500.cov <- apply(na.omit(listit1.1500$cov.treat[1:sims,]), 2, mean)

## average CI length

nls1.500.ci.length <- mean(abs(nls1.500$ub.treat[1:sims,2] - nls1.500$lb.treat[1:sims,2]))
nls1.1000.ci.length <- mean(abs(nls1.1000$ub.treat[1:sims,2] - nls1.1000$lb.treat[1:sims,2]))
nls1.1500.ci.length <- mean(abs(nls1.1500$ub.treat[1:sims,2] - nls1.1500$lb.treat[1:sims,2]))

mle1.500.ci.length <- mean(abs(mle1.500$ub.treat[1:sims,2] - mle1.500$lb.treat[1:sims,2]))
mle1.1000.ci.length <- mean(abs(mle1.1000$ub.treat[1:sims,2] - mle1.1000$lb.treat[1:sims,2]))
mle1.1500.ci.length <- mean(abs(mle1.1500$ub.treat[1:sims,2] - mle1.1500$lb.treat[1:sims,2]))

listit1.500.ci.length <- mean(abs(na.omit(listit1.500$ub.treat[1:sims,2]) - na.omit(listit1.500$lb.treat[1:sims,2])))
listit1.1000.ci.length <- mean(abs(na.omit(listit1.1000$ub.treat[1:sims,2]) - na.omit(listit1.1000$lb.treat[1:sims,2])))
listit1.1500.ci.length <- mean(abs(na.omit(listit1.1500$ub.treat[1:sims,2]) - na.omit(listit1.1500$lb.treat[1:sims,2])))


# SIMULATION 1

## bias

nls2.500.bias <- apply(nls2.500$est.treat[1:sims,], 2, mean) - c(0,1)
nls2.1000.bias <- apply(nls2.1000$est.treat[1:sims,], 2, mean) - c(0,1)
nls2.1500.bias <- apply(nls2.1500$est.treat[1:sims,], 2, mean) - c(0,1)

mle2.500.bias <- apply(mle2.500$est.treat[1:sims,], 2, mean) - c(0,1)
mle2.1000.bias <- apply(mle2.1000$est.treat[1:sims,], 2, mean) - c(0,1)
mle2.1500.bias <- apply(mle2.1500$est.treat[1:sims,], 2, mean) - c(0,1)

listit2.500.bias <- apply(listit2.500$est.treat[1:sims,], 2, mean) - c(0,1)
listit2.1000.bias <- apply(listit2.1000$est.treat[1:sims,], 2, mean) - c(0,1)
listit2.1500.bias <- apply(listit2.1500$est.treat[1:sims,], 2, mean) - c(0,1)

## rmse

nls2.500.rmse <- sqrt(nls2.500.bias^2 + apply(nls2.500$est.treat[1:sims,], 2, var))
nls2.1000.rmse <- sqrt(nls2.1000.bias^2 + apply(nls2.1000$est.treat[1:sims,], 2, var))
nls2.1500.rmse <- sqrt(nls2.1500.bias^2 + apply(nls2.1500$est.treat[1:sims,], 2, var))

mle2.500.rmse <- sqrt(mle2.500.bias^2 + apply(mle2.500$est.treat[1:sims,], 2, var))
mle2.1000.rmse <- sqrt(mle2.1000.bias^2 + apply(mle2.1000$est.treat[1:sims,], 2, var))
mle2.1500.rmse <- sqrt(mle2.1500.bias^2 + apply(mle2.1500$est.treat[1:sims,], 2, var))

listit2.500.rmse <- sqrt(listit2.500.bias^2 + apply(listit2.500$est.treat[1:sims,], 2, var))
listit2.1000.rmse <- sqrt(listit2.1000.bias^2 + apply(listit2.1000$est.treat[1:sims,], 2, var))
listit2.1500.rmse <- sqrt(listit2.1500.bias^2 + apply(listit2.1500$est.treat[1:sims,], 2, var))


## coverage

nls2.500.cov <- apply(nls2.500$cov.treat[1:sims,], 2, mean)
nls2.1000.cov <- apply(nls2.1000$cov.treat[1:sims,], 2, mean)
nls2.1500.cov <- apply(nls2.1500$cov.treat[1:sims,], 2, mean)

mle2.500.cov <- apply(mle2.500$cov.treat[1:sims,], 2, mean)
mle2.1000.cov <- apply(mle2.1000$cov.treat[1:sims,], 2, mean)
mle2.1500.cov <- apply(mle2.1500$cov.treat[1:sims,], 2, mean)

listit2.500.cov <- apply(na.omit(listit2.500$cov.treat[1:sims,]), 2, mean)
listit2.1000.cov <- apply(na.omit(listit2.1000$cov.treat[1:sims,]), 2, mean)
listit2.1500.cov <- apply(na.omit(listit2.1500$cov.treat[1:sims,]), 2, mean)

## average CI length

nls2.500.ci.length <- mean(abs(nls2.500$ub.treat[1:sims,2] - nls2.500$lb.treat[1:sims,2]))
nls2.1000.ci.length <- mean(abs(nls2.1000$ub.treat[1:sims,2] - nls2.1000$lb.treat[1:sims,2]))
nls2.1500.ci.length <- mean(abs(nls2.1500$ub.treat[1:sims,2] - nls2.1500$lb.treat[1:sims,2]))

mle2.500.ci.length <- mean(abs(mle2.500$ub.treat[1:sims,2] - mle2.500$lb.treat[1:sims,2]))
mle2.1000.ci.length <- mean(abs(mle2.1000$ub.treat[1:sims,2] - mle2.1000$lb.treat[1:sims,2]))
mle2.1500.ci.length <- mean(abs(mle2.1500$ub.treat[1:sims,2] - mle2.1500$lb.treat[1:sims,2]))

listit2.500.ci.length <- mean(abs(na.omit(listit2.500$ub.treat[1:sims,2]) - na.omit(listit2.500$lb.treat[1:sims,2])))
listit2.1000.ci.length <- mean(abs(na.omit(listit2.1000$ub.treat[1:sims,2]) - na.omit(listit2.1000$lb.treat[1:sims,2])))
listit2.1500.ci.length <- mean(abs(na.omit(listit2.1500$ub.treat[1:sims,2]) - na.omit(listit2.1500$lb.treat[1:sims,2])))


# SIMULATION 3

## bias

nls3.500.bias <- apply(nls3.500$est.treat[1:sims,], 2, mean) - c(0,1)
nls3.1000.bias <- apply(nls3.1000$est.treat[1:sims,], 2, mean) - c(0,1)
nls3.1500.bias <- apply(nls3.1500$est.treat[1:sims,], 2, mean) - c(0,1)

mle3.500.bias <- apply(mle3.500$est.treat[1:sims,], 2, mean) - c(0,1)
mle3.1000.bias <- apply(mle3.1000$est.treat[1:sims,], 2, mean) - c(0,1)
mle3.1500.bias <- apply(mle3.1500$est.treat[1:sims,], 2, mean) - c(0,1)

listit3.500.bias <- apply(listit3.500$est.treat[1:sims,], 2, mean) - c(0,1)
listit3.1000.bias <- apply(listit3.1000$est.treat[1:sims,], 2, mean) - c(0,1)
listit3.1500.bias <- apply(listit3.1500$est.treat[1:sims,], 2, mean) - c(0,1)

## rmse

nls3.500.rmse <- sqrt(nls3.500.bias^2 + apply(nls3.500$est.treat[1:sims,], 2, var))
nls3.1000.rmse <- sqrt(nls3.1000.bias^2 + apply(nls3.1000$est.treat[1:sims,], 2, var))
nls3.1500.rmse <- sqrt(nls3.1500.bias^2 + apply(nls3.1500$est.treat[1:sims,], 2, var))

mle3.500.rmse <- sqrt(mle3.500.bias^2 + apply(mle3.500$est.treat[1:sims,], 2, var))
mle3.1000.rmse <- sqrt(mle3.1000.bias^2 + apply(mle3.1000$est.treat[1:sims,], 2, var))
mle3.1500.rmse <- sqrt(mle3.1500.bias^2 + apply(mle3.1500$est.treat[1:sims,], 2, var))

listit3.500.rmse <- sqrt(listit3.500.bias^2 + apply(listit3.500$est.treat[1:sims,], 2, var))
listit3.1000.rmse <- sqrt(listit3.1000.bias^2 + apply(listit3.1000$est.treat[1:sims,], 2, var))
listit3.1500.rmse <- sqrt(listit3.1500.bias^2 + apply(listit3.1500$est.treat[1:sims,], 2, var))


## coverage

nls3.500.cov <- apply(nls3.500$cov.treat[1:sims,], 2, mean)
nls3.1000.cov <- apply(nls3.1000$cov.treat[1:sims,], 2, mean)
nls3.1500.cov <- apply(nls3.1500$cov.treat[1:sims,], 2, mean)

mle3.500.cov <- apply(mle3.500$cov.treat[1:sims,], 2, mean)
mle3.1000.cov <- apply(mle3.1000$cov.treat[1:sims,], 2, mean)
mle3.1500.cov <- apply(mle3.1500$cov.treat[1:sims,], 2, mean)

listit3.500.cov <- apply(na.omit(listit3.500$cov.treat[1:sims,]), 2, mean)
listit3.1000.cov <- apply(na.omit(listit3.1000$cov.treat[1:sims,]), 2, mean)
listit3.1500.cov <- apply(na.omit(listit3.1500$cov.treat[1:sims,]), 2, mean)

## average CI length

nls3.500.ci.length <- mean(abs(nls3.500$ub.treat[1:sims,2] - nls3.500$lb.treat[1:sims,2]))
nls3.1000.ci.length <- mean(abs(nls3.1000$ub.treat[1:sims,2] - nls3.1000$lb.treat[1:sims,2]))
nls3.1500.ci.length <- mean(abs(nls3.1500$ub.treat[1:sims,2] - nls3.1500$lb.treat[1:sims,2]))

mle3.500.ci.length <- mean(abs(mle3.500$ub.treat[1:sims,2] - mle3.500$lb.treat[1:sims,2]))
mle3.1000.ci.length <- mean(abs(mle3.1000$ub.treat[1:sims,2] - mle3.1000$lb.treat[1:sims,2]))
mle3.1500.ci.length <- mean(abs(mle3.1500$ub.treat[1:sims,2] - mle3.1500$lb.treat[1:sims,2]))

listit3.500.ci.length <- mean(abs(na.omit(listit3.500$ub.treat[1:sims,2]) - na.omit(listit3.500$lb.treat[1:sims,2])))
listit3.1000.ci.length <- mean(abs(na.omit(listit3.1000$ub.treat[1:sims,2]) - na.omit(listit3.1000$lb.treat[1:sims,2])))
listit3.1500.ci.length <- mean(abs(na.omit(listit3.1500$ub.treat[1:sims,2]) - na.omit(listit3.1500$lb.treat[1:sims,2])))

# SIMULATION 4

## bias

nls4.500.bias <- apply(nls4.500$est.treat[1:sims,], 2, mean) - c(0,1)
nls4.1000.bias <- apply(nls4.1000$est.treat[1:sims,], 2, mean) - c(0,1)
nls4.1500.bias <- apply(nls4.1500$est.treat[1:sims,], 2, mean) - c(0,1)

mle4.500.bias <- apply(mle4.500$est.treat[1:sims,], 2, mean) - c(0,1)
mle4.1000.bias <- apply(mle4.1000$est.treat[1:sims,], 2, mean) - c(0,1)
mle4.1500.bias <- apply(mle4.1500$est.treat[1:sims,], 2, mean) - c(0,1)

listit4.500.bias <- apply(na.omit(listit4.500$est.treat[1:sims,]), 2, mean) - c(0,1)
listit4.1000.bias <- apply(na.omit(listit4.1000$est.treat[1:sims,]), 2, mean) - c(0,1)
listit4.1500.bias <- apply(na.omit(listit4.1500$est.treat[1:sims,]), 2, mean) - c(0,1)

## rmse

nls4.500.rmse <- sqrt(nls4.500.bias^2 + apply(nls4.500$est.treat[1:sims,], 2, var))
nls4.1000.rmse <- sqrt(nls4.1000.bias^2 + apply(nls4.1000$est.treat[1:sims,], 2, var))
nls4.1500.rmse <- sqrt(nls4.1500.bias^2 + apply(nls4.1500$est.treat[1:sims,], 2, var))

mle4.500.rmse <- sqrt(mle4.500.bias^2 + apply(mle4.500$est.treat[1:sims,], 2, var))
mle4.1000.rmse <- sqrt(mle4.1000.bias^2 + apply(mle4.1000$est.treat[1:sims,], 2, var))
mle4.1500.rmse <- sqrt(mle4.1500.bias^2 + apply(mle4.1500$est.treat[1:sims,], 2, var))

listit4.500.rmse <- sqrt(listit4.500.bias^2 + apply(listit4.500$est.treat[1:sims,], 2, var))
listit4.1000.rmse <- sqrt(listit4.1000.bias^2 + apply(listit4.1000$est.treat[1:sims,], 2, var))
listit4.1500.rmse <- sqrt(listit4.1500.bias^2 + apply(listit4.1500$est.treat[1:sims,], 2, var))

## coverage

nls4.500.cov <- apply(nls4.500$cov.treat[1:sims,], 2, mean)
nls4.1000.cov <- apply(nls4.1000$cov.treat[1:sims,], 2, mean)
nls4.1500.cov <- apply(nls4.1500$cov.treat[1:sims,], 2, mean)

mle4.500.cov <- apply(mle4.500$cov.treat[1:sims,], 2, mean)
mle4.1000.cov <- apply(mle4.1000$cov.treat[1:sims,], 2, mean)
mle4.1500.cov <- apply(mle4.1500$cov.treat[1:sims,], 2, mean)

listit4.500.cov <- apply(na.omit(listit4.500$cov.treat[1:sims,]), 2, mean)
listit4.1000.cov <- apply(na.omit(listit4.1000$cov.treat[1:sims,]), 2, mean)
listit4.1500.cov <- apply(na.omit(listit4.1500$cov.treat[1:sims,]), 2, mean)

## average CI length

nls4.500.ci.length <- mean(abs(nls4.500$ub.treat[1:sims,2] - nls4.500$lb.treat[1:sims,2]))
nls4.1000.ci.length <- mean(abs(nls4.1000$ub.treat[1:sims,2] - nls4.1000$lb.treat[1:sims,2]))
nls4.1500.ci.length <- mean(abs(nls4.1500$ub.treat[1:sims,2] - nls4.1500$lb.treat[1:sims,2]))

mle4.500.ci.length <- mean(abs(na.omit(mle4.500$ub.treat[1:sims,2]) - na.omit(mle4.500$lb.treat[1:sims,2])))
mle4.1000.ci.length <- mean(abs(mle4.1000$ub.treat[1:sims,2] - mle4.1000$lb.treat[1:sims,2]))
mle4.1500.ci.length <- mean(abs(mle4.1500$ub.treat[1:sims,2] - mle4.1500$lb.treat[1:sims,2]))

listit4.500.ci.length <- mean(abs(na.omit(listit4.500$ub.treat[1:sims,2]) - na.omit(listit4.500$lb.treat[1:sims,2])))
listit4.1000.ci.length <- mean(abs(na.omit(listit4.1000$ub.treat[1:sims,2]) - na.omit(listit4.1000$lb.treat[1:sims,2])))
listit4.1500.ci.length <- mean(abs(na.omit(listit4.1500$ub.treat[1:sims,2]) - na.omit(listit4.1500$lb.treat[1:sims,2])))


ssize <- c(500, 1000, 1500)
pdf("Figure4.pdf", width = 7, height = 7)
par(mfcol = c(3, 4), mar=c(1.4,1.4,1.4,1.4), oma=c(3,3.5,3,0))

## bias
plot(ssize, c(nls1.500.bias[2], nls1.1000.bias[2], nls1.1500.bias[2]), type = "b", xlim = c(450,1550), ylim = c(0,.26),
     xlab = "Sample Size", ylab = "", pch = 16, lab = c(3, 5, 5))
lines(ssize, c(mle1.500.bias[2], mle1.1000.bias[2], mle1.1500.bias[2]), type = "b")
lines(ssize, c(listit1.500.bias[2], listit1.1000.bias[2], listit1.1500.bias[2]), type = "b", pch = 5, lty = "dashed")
abline(h = 0, lty = "dashed")
legend(x = 900, y = .27, c("LISTIT", "NLS", "MLE"), lty = c("dashed", "solid", "solid"), pch = c(5, 16, 1), bty = "n")

plot(ssize, c(nls1.500.rmse[2], nls1.1000.rmse[2], nls1.1500.rmse[2]), type = "b", xlim = c(450,1550),
     xlab = "Sample Size", ylab = "", pch = 16, ylim = c(0,6), lab = c(3, 3, 5))
lines(ssize, c(mle1.500.rmse[2], mle1.1000.rmse[2], mle1.1500.rmse[2]), type = "b")
lines(ssize, c(listit1.500.rmse[2], listit1.1000.rmse[2], listit1.1500.rmse[2]), type = "b", pch = 5, lty = "dashed")
abline(h = 0, lty = "dashed")
legend(x = 900, y = 6.19, c("LISTIT", "NLS", "MLE"), lty = c("dashed", "solid", "solid"), pch = c(5, 16, 1), bty = "n")

plot(ssize, c(nls1.500.cov[2], nls1.1000.cov[2], nls1.1500.cov[2]), type = "b", xlim = c(450,1550),
     xlab = "Sample Size", ylab = "", pch = 16, ylim = c(.89,1), lab = c(3, 3, 5))
lines(ssize, c(mle1.500.cov[2], mle1.1000.cov[2], mle1.1500.cov[2]), type = "b")
lines(ssize, c(listit1.500.cov[2], listit1.1000.cov[2], listit1.1500.cov[2]), type = "b", pch = 5, lty = "dashed")
abline(h = 0.9, lty = "dashed")
legend(x = 900, y = 1.004, c("LISTIT", "NLS", "MLE"), lty = c("dashed", "solid", "solid"), pch = c(5, 16, 1), bty = "n")

# sim 2
plot(ssize, c(nls2.500.bias[2], nls2.1000.bias[2], nls2.1500.bias[2]), type = "b", xlim = c(450,1550), ylim = c(0,.26),
     xlab = "Sample Size", ylab = "", pch = 16, lab = c(3, 3, 5))
lines(ssize, c(mle2.500.bias[2], mle2.1000.bias[2], mle2.1500.bias[2]), type = "b")
lines(ssize, c(listit2.500.bias[2], listit2.1000.bias[2], listit2.1500.bias[2]), type = "b", pch = 5, lty = "dashed")
abline(h = 0, lty = "dashed")

plot(ssize, c(nls2.500.rmse[2], nls2.1000.rmse[2], nls2.1500.rmse[2]), type = "b", xlim = c(450,1550),
     xlab = "Sample Size", ylab = "", pch = 16, ylim = c(0,6), lab = c(3, 3, 5))
lines(ssize, c(mle2.500.rmse[2], mle2.1000.rmse[2], mle2.1500.rmse[2]), type = "b")
lines(ssize, c(listit2.500.rmse[2], listit2.1000.rmse[2], listit2.1500.rmse[2]), type = "b", pch = 5, lty = "dashed")
abline(h = 0, lty = "dashed")

plot(ssize, c(nls2.500.cov[2], nls2.1000.cov[2], nls2.1500.cov[2]), type = "b", xlim = c(450,1550),
     xlab = "Sample Size", ylab = "", pch = 16, ylim = c(.89,1), lab = c(3, 3, 5))
lines(ssize, c(mle2.500.cov[2], mle2.1000.cov[2], mle2.1500.cov[2]), type = "b")
lines(ssize, c(listit2.500.cov[2], listit2.1000.cov[2], listit2.1500.cov[2]), type = "b", pch = 5, lty = "dashed")
abline(h = 0.9, lty = "dashed")

# sim 3
plot(ssize, c(nls3.500.bias[2], nls3.1000.bias[2], nls3.1500.bias[2]), type = "b", xlim = c(450,1550), ylim = c(0,.26),
     xlab = "Sample Size", ylab = "", pch = 16, lab = c(3, 3, 5))
lines(ssize, c(mle3.500.bias[2], mle3.1000.bias[2], mle3.1500.bias[2]), type = "b")
lines(ssize, c(listit3.500.bias[2], listit3.1000.bias[2], listit3.1500.bias[2]), type = "b", pch = 5, lty = "dashed")
abline(h = 0, lty = "dashed")

plot(ssize, c(nls3.500.rmse[2], nls3.1000.rmse[2], nls3.1500.rmse[2]), type = "b", xlim = c(450,1550),
     xlab = "Sample Size", ylab = "", pch = 16, ylim = c(0,6), lab = c(3, 3, 5))
lines(ssize, c(mle3.500.rmse[2], mle3.1000.rmse[2], mle3.1500.rmse[2]), type = "b")
lines(ssize, c(listit3.500.rmse[2], listit3.1000.rmse[2], listit3.1500.rmse[2]), type = "b", pch = 5, lty = "dashed")
abline(h = 0, lty = "dashed")

plot(ssize, c(nls3.500.cov[2], nls3.1000.cov[2], nls3.1500.cov[2]), type = "b", xlim = c(450,1550),
     xlab = "Sample Size", ylab = "", pch = 16, ylim = c(.89,1), lab = c(3, 3, 5))
lines(ssize, c(mle3.500.cov[2], mle3.1000.cov[2], mle3.1500.cov[2]), type = "b")
lines(ssize, c(listit3.500.cov[2], listit3.1000.cov[2], listit3.1500.cov[2]), type = "b", pch = 5, lty = "dashed")
abline(h = 0.9, lty = "dashed")

# sim 4
plot(ssize, c(nls4.500.bias[2], nls4.1000.bias[2], nls4.1500.bias[2]), type = "b", xlim = c(450,1550), ylim = c(0,.26),
     xlab = "Sample Size", ylab = "", pch = 16, lab = c(3, 3, 5))
lines(ssize, c(mle4.500.bias[2], mle4.1000.bias[2], mle4.1500.bias[2]), type = "b")
lines(ssize, c(listit4.500.bias[2], listit4.1000.bias[2], listit4.1500.bias[2]), type = "b", pch = 5, lty = "dashed")
abline(h = 0, lty = "dashed")

plot(ssize, c(nls4.500.rmse[2], nls4.1000.rmse[2], nls4.1500.rmse[2]), type = "b", xlim = c(450,1550),
     xlab = "Sample Size", ylab = "", pch = 16, ylim = c(0,6), lab = c(3, 3, 5))
lines(ssize, c(mle4.500.rmse[2], mle4.1000.rmse[2], mle4.1500.rmse[2]), type = "b")
lines(ssize, c(listit4.500.rmse[2], listit4.1000.rmse[2], listit4.1500.rmse[2]), type = "b", pch = 5, lty = "dashed")
abline(h = 0, lty = "dashed")

plot(ssize, c(nls4.500.cov[2], nls4.1000.cov[2], nls4.1500.cov[2]), type = "b", xlim = c(450,1550),
     xlab = "Sample Size", ylab = "", pch = 16, ylim = c(.89,1), lab = c(3, 3, 5))
lines(ssize, c(mle4.500.cov[2], mle4.1000.cov[2], mle4.1500.cov[2]), type = "b")
lines(ssize, c(listit4.500.cov[2], listit4.1000.cov[2], listit4.1500.cov[2]), type = "b", pch = 5, lty = "dashed")
abline(h = 0.9, lty = "dashed")

## texts in outer margins
mtext("3 non-sensitive items", side = 3, outer = TRUE, at = 0.23, cex = .9, padj = -1)
mtext("equal probabilities", side = 3, outer = TRUE, at = 0.11, cex = .9, padj = 0.7)
mtext("unequal probabilities", side = 3, outer = TRUE, at = 0.36, cex = .9, padj  = 0.7)

mtext("4 non-sensitive items", side = 3, outer = TRUE, at = 0.75, cex = .9, padj  = -1)
mtext("symmetric probabilities", side = 3, outer = TRUE, at = 0.62, cex = .9, padj = 0.7)
mtext("skewed probabilities", side = 3, outer = TRUE, at = 0.865, cex = .9, padj = 0.7)

mtext("Bias", side = 2, outer = TRUE, at = 0.845, cex = .8, padj = -1.2)
mtext("Root Mean\nSquared Error", side = 2, outer = TRUE, at = 0.5, cex = .8, padj = -.4)
mtext("Coverage of 90%\nConfidence Intervals", side = 2, outer = TRUE, at = 0.15, cex = .8, padj = -.4)

mtext("Sample Size", side = 1, outer = TRUE, line = 0.5, at = 0.125, cex = 0.75, padj = 1)
mtext("Sample Size", side = 1, outer = TRUE, line = 0.5, at = 0.37, cex = 0.75, padj = 1)
mtext("Sample Size", side = 1, outer = TRUE, line = 0.5, at = 0.63, cex = 0.75, padj = 1)
mtext("Sample Size", side = 1, outer = TRUE, line = 0.5, at = 0.88, cex = 0.75, padj = 1)

dev.off()


## #####
## Simulations of the power of the proposed tests to detect design effects

  ## ####
  ## Generate the simulations

## This code was run on a loop, where the number of "yeses to 3 non-sensitive items",
## sample size, and average design effect were varied

library(list)

## number of simulations
sims <- 11000
## number of non-sensitive items
J <- 4
## sample size
n <- 500
## sensitivity parameter
lambda <- c(0, -0.05, 0.05, -0.1, 0.1, -0.2, 0.2, -0.15, 0.15)
## prob for non-sensitive items in the control group
p0 <- c(.25, .5, .75)
## prob for sensitive item
p <- c(.1, .25, .5, .75, .9)

## set alpha for critical value
alpha <- 0.1

tests <- tests.gms <- list()
for(k in 1:length(p)) {
  tests[[k]] <- tests.gms[[k]] <- list()
  for (h in 1:length(p0)) {
    tests[[k]][[h]] <- tests.gms[[k]][[h]] <- matrix(NA, nrow = length(lambda), ncol = sims)
    for (i in 1:length(lambda)) {
      ## prob for non-sensitive items in the treatment group
      p1 <- p0[h] + lambda[i]
      for (j in 1:sims) {
        
        data <- data.frame(y = c(rbinom(n/2, prob = p1, size = J) + rbinom(n/2, prob = p[k], size = 1),
                             rbinom(n/2, prob = p0[h], size = J)), t = c(rep(1, n/2), rep(0, n/2)))
        try(tests[[k]][[h]][i,j] <- ict.test(data$y, data$t, J=J, alpha = alpha, gms = FALSE))
        try(tests.gms[[k]][[h]][i,j] <- ict.test(data$y, data$t, J=J, alpha = alpha, gms = TRUE))
        print(j)
        
      }
    }
  }
}

save(tests, tests.gms, file = paste("TestsN", n,"-J",J,".Rdata", sep = ""))


  ## ####
  ## Figure 5 - "Statistical Power of the Proposed Test to Detect Design Effects"
  ## pg 34

J <- 3

load(paste("TestsN500-J",J,".RData",sep=""))
tests.500 <- tests.gms
tests.std.500 <- tests
load(paste("TestsN1000-J",J,".RData",sep=""))
tests.1000 <- tests.gms
tests.std.1000 <- tests
load(paste("TestsN2000-J",J,".RData",sep=""))
tests.2000 <- tests.gms
tests.std.2000 <- tests

lambda <- c(0, -0.05, 0.05, -0.1, 0.1, -0.2, 0.2, -0.15, 0.15)
## prob for non-sensitive items in the control group
p0 <- c(.25, .5, .75)
## prob for sensitive item
p <- c(.1, .25, .5, .75, .9)

## set alpha for critical value
alpha <- 0.1

##
## simulations start here
##

for(k in 1:length(p)) {
  for (h in 1:length(p0)) {
    for (i in 1:length(lambda)) {
      tests.500[[k]][[h]][i,] <- tests.500[[k]][[h]][i,] <= alpha
      tests.1000[[k]][[h]][i,] <- tests.1000[[k]][[h]][i,] <= alpha
      tests.2000[[k]][[h]][i,] <- tests.2000[[k]][[h]][i,] <= alpha
      tests.std.500[[k]][[h]][i,] <- tests.std.500[[k]][[h]][i,] <= alpha
      tests.std.1000[[k]][[h]][i,] <- tests.std.1000[[k]][[h]][i,] <= alpha
      tests.std.2000[[k]][[h]][i,] <- tests.std.2000[[k]][[h]][i,] <= alpha
    }
  }
}

lambda <- lambda * J

pdf("Figure5.pdf", paper = "special", width = 5, height = 7)

par(mfcol = c(5, 3), mar=c(1.4,1.4,1.4,1.4), oma=c(3,5,2,0), cex = .55, cex.axis = .8)

for(h in 1:3){
  for(k in 1:5){
    power.500 <- power.1000 <- power.2000 <- rep(NA, length(lambda))
    for(i in 1:length(lambda)){
      power.500[i] <- mean(na.omit(tests.500[[k]][[h]][i,])[1:10000])
      power.1000[i] <- mean(na.omit(tests.1000[[k]][[h]][i,])[1:10000])
      power.2000[i] <- mean(na.omit(tests.2000[[k]][[h]][i,])[1:10000])
    }
    
    plot(lambda[order(lambda)], power.500[order(lambda)], xlim = c(-.2*J, .2*J), ylim = c(0,1.01), xlab = "", ylab = "Design", type = "b", pch = 16, axes = F)

    lines(lambda[order(lambda)], power.1000[order(lambda)], type = "b", lty = "dashed")
    lines(lambda[order(lambda)], power.2000[order(lambda)], type = "b", lty = "dotted", pch = 5)

    axis(1, at = c(-.6, -.3, 0, .3, .6))
    axis(2, at = c(0, .25, .5 , .75, 1), labels = c("0.0", "0.25", "0.5", "0.75", "1.0"))

    box(which = "plot")

    if(k==1) legend(x = -.05, y = 1.05, c("n=500", "n=1000", "n=2000"), lty = c("solid", "dashed", "dotted"), pch = c(16, 1, 5), bty = "n", cex = .9)
    
  }
}
  
## texts in outer margins
mtext(expression(paste("Expected Number of Yeses to ", 3, " Non-sensitive Items ",E(Y[i](0)))), side = 3, outer = TRUE, at = 0.51, cex = .7, padj = -0.5)
mtext("0.75", side = 3, outer = TRUE, at = 0.17, cex = .7, padj = 1)
mtext("1.5", side = 3, outer = TRUE, at = 0.5, cex = .7, padj  = 1)
mtext("2.25", side = 3, outer = TRUE, at = 0.84, cex = .7, padj  = 1)

mtext(expression(paste("Probability of Answering Yes to the Sensitive Item Pr(Z", scriptstyle("i, J+1"), ")")), side = 2, outer = TRUE, at = 0.5, cex = .7, padj = -3.5)
mtext("*", side = 2, outer = TRUE, at = 0.695, cex = .7, padj = -5)
mtext("0.1", side = 2, outer = TRUE, at = 0.9, cex = .7, padj = -2.4)
mtext("0.25", side = 2, outer = TRUE, at = 0.7, cex = .7, padj = -2.4)
mtext("0.5", side = 2, outer = TRUE, at = 0.5, cex = .7, padj = -2.4)
mtext("0.75", side = 2, outer = TRUE, at = 0.3, cex = .7, padj = -2.4)
mtext("0.9", side = 2, outer = TRUE, at = 0.1, cex = .7, padj = -2.4)

mtext("Statistical Power", side = 2, outer = TRUE, at = 0.9, cex = .6, padj = -0.9)
mtext("Statistical Power", side = 2, outer = TRUE, at = 0.7, cex = .6, padj = -0.9)
mtext("Statistical Power", side = 2, outer = TRUE, at = 0.5, cex = .6, padj = -0.9)
mtext("Statistical Power", side = 2, outer = TRUE, at = 0.3, cex = .6, padj = -0.9)
mtext("Statistical Power", side = 2, outer = TRUE, at = 0.1, cex = .6, padj = -0.9)


mtext(expression(paste("Average Design Effect ",Delta)), side = 1, outer = TRUE, at = 0.17, cex = .6, padj = 1.2)
mtext(expression(paste("Average Design Effect ",Delta)), side = 1, outer = TRUE, at = 0.5, cex = .6, padj  = 1.2)
mtext(expression(paste("Average Design Effect ",Delta)), side = 1, outer = TRUE, at = 0.84, cex = .6, padj  = 1.2)

dev.off()

